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Abstract 

An algorithm for the generation of non-uniform unstructured grids on ellipsoidal geometries is described. 
This technique is designed to generate high quality triangular and polygonal meshes appropriate for general 
circulation modelling on the sphere, including applications to atmospheric and ocean simulation, and nu- 
merical weather predication. Using a recently developed Frontal-Delaunay-refinement technique, a method 
for the construction of high-quality unstructured ellipsoidal Delaunay triangulations is introduced. A dual 
polygonal grid, derived from the associated Voronoi diagram, is also optionally generated as a by-product. 
Compared to existing techniques, it is shown that the Frontal-Delaunay approach typically produces grids 
with near-optimal element quality and smooth grading characteristics, while imposing relatively low com- 
putational expense. Initial results are presented for a selection of uniform and non-uniform ellipsoidal 
grids appropriate for large-scale geophysical applications. The use of user-defined mesh-sizing functions to 
generate smoothly graded, non-uniform grids is discussed. 
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1. Introduction 

Planetary climate simulation and numerical weather predication are two of the most computationally 
challenging tasks currently under consideration. General circulation models typically involve the solution of 
a large coupled set of non-linear transport equations on very high-resolution spatial grids. The complexity 
of these models is often significant, with the discrete systems typically involving dozens of independent 
variables, including fluid momentum, thermodynamic state, such as temperature, salinity and precipitable 
water, atmospheric aerosols and oceanic biogeochemical tracers. While many such models have previously 
been developed using structured discretisations of the underlying spherical geometry, recent interest in the 
use of unstructured general circulation models necessitates the development of effective unstructured mesh 
generators that cater to this class of numerical scheme. This study investigates the applicability of a recently 
developed Frontal-Delaunay surface meshing algorithm for this task. 

1.1. Related Work 

While a simple structured mesh for the sphere can be obtained by building a uniform discretisation in 
spherical coordinates, the resulting ‘lat-lon’ grid is inappropriate for numerical simulation, incorporating 
a pair of grid singularities at the poles. A majority of current generation general circulation models, (i.e. 
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Figure 1: Semi-structured meshing for the sphere, showing a cubed-sphere grid on the left, and an icosahedral grid on the 
right. Both grids are generated with equivalent mean edge lengths. 



M), are instead based on a semi-structured discretisation known as the cubed-sphere. According to this 
methodology, the spherical surface is decomposed topologically into an array of six quadrilateral ‘faces’, 
with each face in turn decomposed into a simple structured curvilinear grid. Each face of the cube is 
projected onto the adjacent spherical geometry. The cubed-sphere topology is illustrated in Figure [I] In 
addition to the cubed-sphere configuration, a second class of semi-structured spherical grid can be obtained 
through icosahedral-type decompositions. In such cases, the primary grid is a triangulation, though the dual 
polygonal grid is often also used. Icosahedral-type grids are also illustrated in Figure [l] 

Recent general circulations models, specifically the Model for Predication Across Scales (MPAS) [31 
110 , has focused on the use of locally refined, unstructured polygonal grids to achieve multi-resolution 
representations of both atmospheric and oceanic dynamics. The current study explores the development of 
an algorithm for the generation of high resolution spheroidal Delaunay triangulations and associated Voronoi 
diagrams appropriate for such climate models. Note that in addition to the high-resolution requirements, 
the sensitivity of the underlying numerical schemes employed in general circulation models necessitates the 
construction of meshes with near-optimal mesh quality. 


2. A Frontal-Delaunay- Refinement Grid- Generator 


The task is to generate very high-quality meshes for planetary atmospheres and/or oceans. Grids are 
a semi-structured ‘2.5D’ configuration, in which a two-dimensional unstructured grid generated over the 
spheroidal surface is extruded in the radial direction to form a semi-structured ‘stack’ of fluid layers. The 
layer-wise grid is set-up to follow the planetary topography - typically represented as a set of elevations 
defined from a reference spheroid. In a general form, this reference surface can be expressed in terms of an 
axis-aligned triaxial ellipsoid: 
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where the Xi s are the Cartesian coordinates in a locally aligned coordinate system, the point (c±, C 2 , C 3 ) is 
the centre of the ellipsoid, and the scalars (?~i , s ^ 3 ) are its principal radii. 
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2.1. Mesh- size Functions 

Local mesh density can be controlled via a user-specified mesh-size function h(pc) defined on a background 
grid. Continuous sizing values are computed using an interpolation scheme. In the current study, mesli-size 
functions are composed using structured ‘lat-lon’ grids and bilinear interpolation providing support for 
mesh-size functions derived from high-resolution satellite altimetry data [6|. 

2.2. Delaunay Mesh Generation 

A high-quality triangular surface mesh is generated on the ellipsoidal reference surface (pi) using a 
‘Frontal-Delaunay’ variant of the conventional restricted Delaunay-refinement algorithm m misi usi m]- 
This technique is described by the author in detail in mm and differs from standard Delaunay-refinement 
approaches in terms of its placement of Steiner vertices. Specifically, the Frontal-Delaunay variant employs 
various off-centre techniques [H}[l5], designed to position refinement points such that element- quality and 
mesh-size constraints are satisfied in a locally optimal fashion. Previous studies have shown that such an 
approach typically leads to substantial improvements in mean element-quality and mesh smoothness. Addi- 
tionally, it has been demonstrated that this method inherits much of the theoretical robustness of standard 
Delaunay-refinement techniques - offering guaranteed convergence, topological correctness, and a minimum 
angle guarantee. 

Given a user-defined mesh-size function h(x) and an upper-bound on the element radius-edge ratios p, 
the Frontal-Delaunay algorithm proceeds to sample the ellipsoidal surface by refining any surface triangle 
that violates either the mesh-size or element-quality constraints. Refinement is accomplished by inserting 
a new Steiner vertex at the off-centre point associated with a given element. Refinement continues until 
all constraints are satisfied. Upon termination, the resulting surface mesh is guaranteed to contain nicely 
shaped surface triangles, satisfying both the radius-edge constraints p(fi) < p and the mesh-size limits 
M/i) < M x /) f° r a ll surface triangles /* in the mesh. 

As a restricted Delaunay-refinement approach, a full three-dimensional Delaunay tetrahedralisation 
Del(X) is incrementally maintained throughout the surface meshing phase, where X £ R 3 is the set of vertices 
positioned on the surface of the ellipsoid geometry. The set of restricted surface facets Del |s(A) that conform 
to the underlying ellipsoidal surface are expressed as a subset of the tetrahedral faces Del |s(Af) C Del(X). 
In an effort to minimise the expense associated with maintaining the full-dimensional topological tessellation, 
an additional scaffolding vertex x s is initially inserted at the centre of the ellipsoid. This has the effect of 
simplifying the resulting topological structure of the mesh, with the resulting tetrahedral elements forming 
a simple configuration in which they emanate radially outward from x s . 

2.3. Mesh Optimisation 

While the Frontal-Delaunay-refinement algorithm typically produces surface triangulations of very high- 
quality, these tessellations can be further improved through a mesh optimisation step. In this study, a 
spring-based relaxation strategy, adapted from the DISTMESH algorithm of Persson and Strang [16], is used 
to optimise vertex positions. Given an initial vertex distribution x°, generated using the Frontal-Delaunay- 
refinement algorithm described previously, the relaxation procedure seeks updates to the vertex positions in 
order to achieve an approximate force- equilibrium between the mesh vertices 

x fc+1 = x fc + At F(x fc ) , Fj(x fe ) = ^ F itj (x k ) ■ e xJ , i^(x fe ) = f (J ii3 -(x fc ) - /^(x fc )). (2) 
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Here IZi is the set of vertices adjacent to the i-tli vertex, lij(x. k ) is the length of the edge between the 
vertices x, ; and x ; , hij(x. k ) is the value of the mesh-size function computed at the associated edge midpoint, 
and e, j is the corresponding unit edge direction vector. The spring constant /3 is set to unity in this study. 
Noting that Q does not constrain vertices to remain on the ellipsoidal reference surface, an additional 
projection operator is emloyed within each iteration to enforce the surface constraints exactly. Incremental 
modifications to the underlying Delaunay mesh topology are enacted throughout the iteration process, 
ensuring that the mesh remains fully Delaunay. 
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Figure 2: A multi-resolution polygonal grid of the world ocean, generated using the Frontal-Delaunay- refinement algorithm 
described previously. The non-uniform mesh-size function is computed based on the local gradient in bathymetry, ensuring 
that regions of rapid change are better resolved. 



Q.f. Construction of Voronoi Cells 

Following the generation and optimisation of the surface triangulation, a dual polygonal grid is derived 
from the associated Voronoi complex. The surface Voronoi diagram contains a grid-cell for each vertex in 
the underlying surface triangulation, where each grid-cell is a closed, non-planar convex polygonal region 
that connects the centres of the Surface-Delaunay-Balls associated with the triangles adjacent to the central 
vertex. 


3. Results &: Discussions 

Initial investigations have been conducted using the Frontal-Delaunay-refinement approach described 
previously. A range of very high-resolution spherical grids have been generated using both uniform and 
non-uniform mesh-size functions. A medium resolution polygonal mesh of the global ocean, generated using 
a bathymetrically refined mesh-size function is shown in Figure [2] Future invetigations will analyse the 
effectiveness of solution- adaptive grid generation schemes, including, for instance, those based on time- 
averaged distributions of Eddy Kinetic Energy (EKE). 

While the performance of the present algorithm appears to be reasonable, generating in excess of several 
million elements per minute using standard desktop computing infrastructure, future work will focus on fur- 
ther improvements to performance, including the investigation of dedicated techniques for the computation 
of Spherical Delaunay Triangulations (SDT’s) |TT) , and the development of parallel meshing techniques. A 
detailed comparison with the existing Parallel Spherical Centroidal Voronoi Tessellation algorithm (MPI 
SCVT) of Jacobsen et al. [THj will also be conducted. 
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